function [Sa,Sb,Sp] = pmvdr_PMCC(x,p)

x = x-mean(x); %avoid DC
Nfft = 400;
%tic
[m,n] = size(x);
N = p;
R = ifft(abs(x).^2);
R = R./m; % Biased autocorrelation estimate

[a,e] = levinson(R,N);
lenx = length(x);
%toc

%tic
for k = 0:p,
    i = [0:p-k];
    mu(k+1) = sum((p+1-k-2*i).*a(i+1).*conj(a(i+1+k)));
end
%toc

mu = [conj(mu(end:-1:2)) mu]/e;
lmu = length(mu);

    Wmat = exp(-sqrt(-1)*2*pi/Nfft*[0:Nfft/2-1]'*[0:(lmu-1)/2]);
    Pwa = Wmat*a';
    Sa = e./abs(Pwa).^2/(p+1); %/lenx;

    Wmat = exp(-sqrt(-1)*2*pi/Nfft*[0:Nfft/2-1]'*[-(lmu-1)/2:(lmu-1)/2]);
    Pwb = real(Wmat*mu');
    Sb = 1./Pwb;

    Sp = abs(fft(x,Nfft))/Nfft;
    Sp = Sp(1:Nfft/2).^2;    
    
end
